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I. INTRODUCTION 



The Hartree-Fock-Bogoliubov (HFB) approximation was outlined more than twenty years ago [1] for infinite systems 
and almost immediately was introduced in nuclear physics [2] . In the case of infinite systems the HFB procedure is 
well studied and the character of the wave-functions (w.f.) is well understood. However, in the case of finite systems 
(nuclei), things are not so clear. The HFB equations in the case of nuclei are meaningful provided the boundary 
conditions for the single quasiparticle (s.qp.) w.fs. are correctly formulated in order to describe a genuine finite state. 
The goal of the present paper is to provide the correct formulation of the HFB approximation in the case of finite 
systems. 

It is well known that pairing correlations always appear whenever a pole (i.e. a bound state) is present in the 
two-body Green function of the many-body system [3]. In this case, corrections to the single particle (s.p.) Green 
function of the type 



- £ 

h 




give rise to diagrams of the type 



p h P 

> O < O > (!) 



when only the contribution of the pole is taken into account [3] . Here 2A stands for the energy of the two-particle 
bound state and e p and Sh are the particle and hole energies respectively. The process represented by diagram (1) 
leads to a mixing between particle and hole states and as a result to a smearing of the Fermi surface [3-5]. This 
mixing has the special feature that a hole (particle) can transform into a particle (hole), due to the presence of the 
pair condensate, provided their energies are related by the energy conservation law 

e p + Eh = 2A. 

Whenever the energy of the hole state is less than 2A the corresponding particle state to which the hole state is 
coupled lies in continuum. This situation formally resembles the case of an electron in a very strong field [6] or the 
case of a bound state embedded in continuum [7] . Consequently, a sufficiently deep hole state becomes unstable with 
respect to decay into a particle state by an interaction with the pair condensate. The same thing happens in the case 
of a deep hole decaying into a less deeply bound hole, with the excitation of a phonon, which can further decay by 
particle emission. This process is described by the diagram 




where the wavy line represents a phonon. 

The physical situation is thus not new. New is the fact that due to the presence of the pair condensate, some 
hole states acquire a width and therefore the s.qp. w.fs. can no longer be treated as corresponding to bound states, 
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as has hitherto been done, but as continuum states. How to introduce this property of the s.qp.w.fs. into the HFB 
approximation for finite systems and therefore how to define the boundary conditions correctly is the main goal of 
the present discussion. 



II. HFB EQUATIONS FOR FINITE SYSTEMS 



In this section the well-known HFB equations will be derived for the sake of completeness. The emphasis will be on 
the correct definition of s.qp.w.fs. so as to describe genuine finite systems, i.e. systems with finite matter distribution. 
The forces between particles will not be specified, except for some general properties, like the finite range. 

By analogy with the usual HF approximation, the HFB ground state w.f. | 0) is defined as the vacuum for the 
fermi quasiparticles [1,2,4,5] 



and tp(x) and tp^x) stand for field operators for annihilation and creation of a particle with space-spin coordinates 
x = (r , er) , which satisfy the usual anticommutation relations 



ai | 0) = , (0 | 0) = 1 



(2) 



where 




(3) 




(4) 



{ai,aj} = {at, at} = 0, 



one easily obtains the relations 




(5) 




(6) 




The constraints (5-6) ensure the unitary character of the transformations (3). 
The total energy of the many-body system and the mean number of particles are 



£ = {0\H\0) 



and 



N = (0 | JV | 0) 



where H and N = J dxtp^ (x)tp(x) stand for the hamiltonian and the number operator in the second quantization 
representation. 

The mean values for the energy and for the particle number can be expressed through the densities 
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p(x,y) = (0 | ^ + ( V mx) | 0) = p*(y,z) = J>(zK(y), (7) 

i 

*(x, y) = (0\ i>{y)i>{x) | 0) = -$( ? y, a) = Y, < x ) u i(y)i ( 8 ) 



in the following way 



£ = Tr(7» + l -Tr{Vpp) a + ±Tr(V***), (9) 

N = Tr(p), (10) 

if only two-body interactions are present. T and F stand for the kinetic energy and the two-body interaction 
respectively. A shorthand notation was used for the traces in Rels. (9) and (10). 

The HFB equations for the two component s.qp.w.fs. {vi(x), u*(x)} are derived from the stationarity condition of 
the total energy (9) under the subsidiary conditions (10) and (5). These equations are 

J h(x,y)vi(y)dy - Xvi(x) + J A(x,y)u*(y)dy = E^x) (11) 
J ' £J{x,y)vi{y)dy - J h*(x,y)u*(y)dy + Xu*(x) = E l u l (x) (12) 



where 



is the chemical potential, and 



and 



5£ n 

X = 777 < 
SN 



h ^ y) = J^y) =k *^ X) 



are the s.p. hamiltonian and the pairing field respectively and Ei stands for the s.qp. energies. 

When performing the summations in Rels. (7) and (8) one must include only those solutions of the nonlinear system 
(11-12) with Ei < 0, which define the operators a\. The solutions with Ei > correspond to operators a\. 

Let us analyse in more detail these equations in the case of a finite system. The problem which arises is the meaning 
of the normalization condition (5); namely, if the right hand side of Rcl. (5) should be a (5-function or a Kronecker 
symbol, as is usually the case [2,5,8,9], in complete analogy with the HF approximation. 

A finite system is characterized by a finite matter distribution and therefore by a finite range of the s.p. field (except 
in the case of the Coulomb interaction). Naturally, being determined by the matter distribution, the pairing field 
must also have a finite range (an infinite range of the pairing field can only occur if the system under consideration is 
unstable with respect to two particle decay.) The problem is determining the mechanism which leads to finite matter 
distribution when s.qp. states no longer have the character of discrete states, as was discussed in the Introduction. 

When pairing is turned off (i.e. A = 0) the spectrum of Eqs. (11-12) looks like the one show in Fig. 1. The left 
hand side corresponds to the spectrum of Eq. (11), while the right hand side corresponds to Eq. (12). 
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FIG. 1. The spectra of Eqs. (11-12) when the pairing field is turned off (A == 0). The left hand side corresponds to the 
spectrum of Eq. (11), while the right hand side corresponds to Eq. (12) . The hatched regions represent the continuum and 
the lines correspond to discrete states. 

Upon turning the pairing on, the two spectra mix and the discrete states with energies outside the interval 

\<E<-\ (13) 

will lie in the continuum. Only the states with energy within the interval (13) will preserve the bound state character. 

The continuum part of the spectrum for E < A and E > —A can disappear only if the s.p. potential is finite and 
positive at infinity as in the case of a harmonic oscillator. 1 There is no physical reason for this to be true in the case 
of finite systems. 

A glance at Fig.l. reveals why pairing correlations lead to a significant increase of the level density in the vicinity 
of the Fermi surface. Pairing appears when the chemical potential has a value within a shell, as illustrated in Fig. 
1. In such a case, the total number of hole (the left hand side of Fig. 1) and particle (the right hand side of Fig. 1) 
states, which is practically equal to the number of s.qp. states with pairing included, is almost double the number of 
s.p. states in HF approximation. 

In order to determine what happens in the case of a bound state with an energy E, when pairing is turned on, we 
shall use perturbation theory. For the sake of simplicity, the nucleus will be assumed to be spherical and the pairing 
field real A = A*. Also, the spin and angular variables are assumed to be already separated from Eqs. (11) and (12) 
and the corresponding geometrical factors included in the definition of the single particle fields. From Eq. (11) one 
easily obtains in the vicinity of a bound state the relation 



1 The discrete character of the spectrum for a limited number of states with E < \(E > —A) can also be preserved if [h, A] = 0, 
a condition which is not fulfilled in fact. This condition is satisfied in the case of constant pairing approximation (A = const.), 
an approximation which leads to an unphysical density distribution. 
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Mr) = (^T^ A ^) W = Mr) ( ^ Al E U Q E) = n^(E)Mr), (14) 

where 

{h-\-E )<f> = , (0o I 0o) = 1. 
The other component of the single quasi particle wave function becomes then 

u E {r) = C(E)u 0E (r) n^ 2 {E) ( ^ + \ _ fe A0 o ) (r), (15) 

where 

(h — X + E)u QE = 0, {u 0E \ u 0E >) = 6(E - E') (16) 
and C(-E') is a normalization constant. From Rcls. (14) and (15) one obtains 

n i/2 {E) = C(E)(<f> a | A | u 0E ) 



E - E + (0 O | A srk r A 



'-£+A-/i 



In order to determine the normalization constant C(E), we shall use the representation of the Green function 
through regular and irregular solutions of Eq. (16) 

G(r, r', —E + X) = ( - - - ) (r,r') = ^e^Xe^) 
V ' \-E + \-h) y ' h 2 W{u 0E , X e) 

where W(uqe,Xe) stands for the Wronskian 



W(u 0E ,Xe) = ~^2- 
nn 

The asymptotic behavior for the regular u 0E and irregular \e solutions is given by 



u 0E (r) = \j Tr - r sin(kr + S 0E ), 



2m 
2m 



xe(t) = \l r cos ( kr + $oe), 



where k 2 — 2m(—E + X)/fi 2 . The asymptotic behavior of the u-component is 



u E (r) ^ J-^L[C(E) sin(fcr + S 0E ) + -Kn 1,2 {E){<j) \ A \ u 0E ) cos(fcr + 5 0E )}, 
n nk 



and therefore 



C 2 (E] = [E - E Q + (0q I A^-gj^A | O )] 2 

[E - E + < 0o I A _ E ^ x _ h A | 0o >} 2 +tt 2 |< 0o i A | u 0E >| 4 ' 1 ' 

Consequently, the normalized solutions are 

v E (r)=n 1 / 2 (E)Mr) (18) 
u E (r) = C(E)u 0E (r) n^ 2 (E) ( _ - -L _ A0 O ) (r) , (19) 

where 
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n(E) 



1 \Y{E) 

ir (E-E -SE) 2 + \T 2 (E) ' 



(20) 



T(E) = 2n\(<f> | A | u 0E )\ 2 = 2tt- 



(</> I [ft, A] | u 0E ) 



(E + S ) 2 



<5£ : 



A- 



1 



-E+X-h 



A^o) 



and 



/ n(£)d£ 

«/ OO 



(21) 



The matrix elements in the above relations can be calculated without any significant loss of accuracy for E = E a . 

As one can observe, due to the coupling with the continuum, the bound state spreads over the entire spectrum (see 
Rcl. (21)). Now the quantity n(E) has thus to be interpreted as the occupation number probability density over a 
unit energy interval. 

The solution is formally equivalent to the solution of a coupled channel problem with a bound state embedded in 
the continuum [7]. It displays a well-defined resonant character with a width T(E ) and a shift 5E. The case of a 
resonant state can be treated in a similar way. 

Far from the resonance energy Eq, the amplitude of the w-component is very small, while the w-component is 
practically equal to the unperturbed solution u e- In the vicinity of the resonance E , the amplitude of the u- 
componcnt increases significantly inside the potential well (see Rels. (18), (20)), while the phase of the w-component 
changes by ir. The phase of the w-component is 



$e = 3qe ~ arctan 



T(E) 



2(E - E - 5E) 



The solution just described is characterized by the fact that the w-component is square integrable, even though the 
single quasi-particle wave function solution (ve(x), u* e {x)) represents a continuum state. The question is: does this 
feature hold true for the selfconsistent solution of the HFB equations? 

The density distribution p satisfies Rel. (10) (i.e. the diagonal part of p is integrable), while the density of the pair 
condensate satisfies the condition 



/ 



*{x,y) | 2 = Tr[p{\ - p)] < N, 



(22) 



which can be easily derived by means of Rels. (5-6) [1,4,5,8]. Therefore, one can expect that p and $ both fall 
down quickly enough outside the system. The density distribution p determines the asymptotic behavior of the single 
particle selfconsistent potential U, while the density $ defines the pairing potential A = V&/2. 

The nonlocality of these potentials is governed by the range of the two-body interaction V, assumed to be finite. 

We shall now show that the following asymptotic behaviors take place 



U <~ p ~ O exp 



(23) 



and 



A ~ $ ~ O \ exp 



/ 2 2m|A| r 



(24) 



when R= (r + r')/2 tends to infinity and s = r — r' remains finite (practically of the order of the range of V). (In 
the above Rels. (23-24) the symbol O means that the quantities on the left hand side behave like the corresponding 
arguments of the O-function.) As one can observe, the pairing field A has a longer tail than the s.p. selfconsistent 
field U. 

Using Rels. (23) and (24) one can show, using HFB Eqs. (11-12), that the v- and w-components behave asymp- 
totically as 
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v(r)^e X p[-J- T \E + X\r\, (25) 



u{r)~exp\-^\E-\\rj, (26) 

if 

A < E < 0, (27) 

and 



V (r)~exp^-^|2A|rJ, (28) 

u(r) ~ 0(1), (29) 

if 

E < X. (30) 

Outside the potential well, for energies in the interval (27), the two Eqs. (11-12) decouple and the corresponding 
asymptotic behavior of the v— and u-components of the single quasi-particle wave function is determined by the 
"energies" A + E and A — E, respectively. For energies in the interval (30,) the asymptotic behavior of the v- 
component is governed by the inhomogeneous part of the Eq. (11) (i.e. by the term Ait*), which cannot be neglected 
in this case, as was possible for energies in the interval (27). On the other hand, the term A^v falls down exponentially 
and does not influence the asymptotic behavior of the ^-component. 

The asymptotic behavior of the w-component is fully determined by the "energy" —E + A, which is positive in the 
interval (30). 

Now, if one takes into account the definitions of densities p and $ (Rcls. (7) and (8) respectively), one can easily 
notice that the asymptotic behavior of the v- and w-components (see Rels. (25-30)) completely agrees with the 
asymptotic behaviors (23) and (24). This means that the corresponding asymptotic behaviors are selfconsistent. It is 
physically natural to expect that the asymptotic behavior of the density p is controlled by the chemical potential A, 
i.e. by the energy of the least bound particle. The density of the pair condensate $ can be interpreted as the wave 
function of a bound state of two interacting particles in an external field with energy 2A. Using HFB Eqs. (11) and 
(12) and the definition (8), one can show that the density $ satisfies the equation [1] 

(h(l) + h{2) + w(l, 2) - 2A)$(1, 2) = (31) 

other terms negligible outside the system . 
From this equation it follows that at large distances $ behaves as 

exp (—ari - j3r 2 ) , 

where 

2m, 

T 2 

which also agrees with Rel. (24). 

Strictly speaking, this asymptote is correct only outside the range of V . If,however, one takes into account the fact 
that a system of two identical nucleons does not have bound states, the asymptote is valid everywhere outside the 
s.p. potential well. 

Summing up, in the normalization conditions (5) for the single quasi-particle wave functions, the right hand side 
has to be interpreted as a Kronecker symbol if A < E < —A and as a Dirac 5-function if | E \>\ A | (as it is well 
known [1-5,8] the system (11-12) has the property that if {vi, u*, Ei} is a solution, then {u i} v* , —Ei} is a solution as 
well). Furthermore, for E < the v-component of the single quasi-particle wave function is always square integrable 
and its norm has to be interpreted as the occupation number probability. On the other hand, the relation 



V+/3 2 ) = 9r 2A > 



1 — fii = 1 — J | Vi(x) \ 2 dx = J | Ui{x) \ 2 dx 



is valid only for A < E < — A. If the energy E is outside this interval, the integral J \vi(x)\ 2 dx should be interpreted 
as an occupation number probability density per unit energy interval. 
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III. SOME SIMPLE EXAMPLES 



This section is devoted to some simple examples of HFB equations, which although somewhat unrealistic, lead to 
a better understanding of various issues arising while solving Eqs. (11-12). 



A. Constant pairing approximation A = constant 

This approximation is also known as the BCS approximation [10]. The s.qp.w.f. in this case is 

v(r) = nJc S <j)(r), 
u(r) = VI - n B cs<P( r )i 

where 

( "- e, * = °' -°-l('- t >.-',y + j < 32) 

and 



E = -y/(e- A) 2 +A 2 . 

Usually, the pairing field A = const is taken to be different from zero in a limited energy interval around the Fermi 
surface. However, there is no recipe to determine this energy band in a unique way. It is obvious that the density p 
cannot be integrable if solutions belonging to the continuum part of the spectrum of the Eq. (32) are included (i.e. 
when the pairing field is nonvanishing for such states). Furthermore, if one considers that A is acting for all energies, 
then the density $ is 

A 



$(x,y) = 4>i(x)4>* (y) , 

V y ' Y K J 2y/(ei- A) 2 + A 2 



A 



^/(si - A) 2 + A 2 e>-\ 



where G(x,y,X) stands for the single particle Green function of the Eq. (32). As is well known, the Green function 
diverges like | r — r' | _1 when r — ► r'. Therefore, a finite density $ cannot be defined in this case because of this 
divergence. One notices that the divergence is not logarithmic, as is usually stated in textbooks [3,10]. A local pairing 
field corresponds to a zero range two-body interaction. Then, as one can easily show by using Eq. (31), the density 
$ will always be singular for coinciding arguments. For the entire HFB procedure to be meaningful, the two-body 
interaction, which is responsible for pairing, must have a finite range. 



B. Square well single particle potential and square well pairing potential 

We assume that 

U(r) = U 6(R - r), U < 0, 
A(r) = A 9(R - r), 

and look for solutions of Eqs. (11-12) with zero orbital momentum I = 0. 
For r < R the solution reads 

Vin{r) = C + sin k + r + C- sin fc_ r, 

Uin (0 = P+ C + sin k + r + [3- CL sin fc_ r, 

where 
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and C± stand for some constants, which have to be determined from matching the interior with the exterior solutions 
and normalization. 
For r > R 

h 2 

v ou t(r) = v exp(-k r), if - — k 2 = —{E + A) > 0, 

Am 

h 2 

u out {r) = uq exp(-ik' r) if t^^q = ~{~E + A) > 



or 



if 



u ou t{r) = Mi sin k' r + u 2 cos k' r 



h 2 

—k' 2 = —E + A > 
2m 

where v ,u ,ui and u 2 have to be determined from matching and normalization. 

If A < E < —A, the spectrum is discrete and the energies have to be determined from the matching condition 

/3+(k- cos k_R + ko sinfc_i?)(fc + cos k + R + k' Q sin k + R) ^ 
(3-(k- cos k^R + k' sin k-R)(k+ cos k + R + k sin k+R) 

If E < A , then one deals with a state lying in continuum. In contrast to the case of constant pairing, one now has 
[h, A] 7^ and sufficiently deeply bound hole states acquire width. 

Inside the potential well, the two components of the s.qp.w.f. form a superposition of two s.p.w.f. with energies 
equal approximately to E + A + Uo and — E + A + Uq, respectively (as a rule A is very small and can be neglected 
in the determination of k±). 

In the vicinity of a hole state (the left hand side of Fig. 1) C+ has a zero and one can show that 



v in(r) ^ a ^ 2E I uses ^ -, 

- P+ = -r- = \ r. , riBCS = 1, 



u in (r) A V 1 ~ n BCS 

similar to the BCS approximation. This relation holds only inside the potential well and it is not valid outside it. Far 
from the resonance this ratio becomes 



Uin(r) ^ a _ 2E I - uses „ , 

tt = P+ = -T- = \ 7 n BCS < 1- 

Vi n {r) A V rises 

Unlike the BCS approximation, the radial behaviors of the v- and u-components are no longer identical. 

In this case, the density p has the correct asymptotic behavior, but due to the local character of the pairing potential, 
the density $ has the same divergence as in the case of BCS approximation. . 

C. Surface pairing A = A 5(r — R) 

This is another approximation which has been used in nuclear physics. The solution of Eqs. (11-12) is now 

v(r)=A G(r,R,E + X)u(R), 

u(r) = Ku (r) - AgG(r, R, -E + X)G(R, R,E + X)u(R), 

where 

KUo(R) 



u{r) 



1 + A 2 G(R, R, -E + X)G(R, R 7 E + X)' 
k is a normalization constant and uq is given by the equation 
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(h - A + E)u Q = 



(u has to be included only for E < A) 

The density p has a correct asymptotic behavior, but the same problems arise with the density <I> as above. The 
density of the pair condensate $ is singular for r = r' (see Eq. (31)). 

Even though the examples discussed here have little in common with the real selfconsistent solution of Eqs. (11-12), 
in our opinion however, they lead to a deeper understanding of the structure of the HFB equations in the case of finite 
systems. Especially instructive in this sense is the role played by the nonlocality of the pairing field and consequently 
by the range of the two-body forces. 



IV. CONCLUSION 



We have examined the HFB approximation in the case of finite systems, when the two-body interaction between 
particles has a finite range. Special attention was paid to the asymptotic behavior of the s.qp.w.fs. It was shown that 
the s.qp. states located in spectrum sufficiently far from the Fermi surface have the character of continuum states. 
E.g. the deep hole states acquire a width corresponding to the decay into a particle state and the pair condensate. 
This width has to be interpreted as a contribution to the imaginary part of the s.p. optical potential. 

Even though most of the s.qp.w.fs. lie in continuum, the matter distribution is nevertheless finite. The same holds 
true for the density of the pair condensate, which is finite as well. 

These features of the general solutions of the HFB equations have to be included in any HFB calculations. To what 
extent the available HFB results [8,9] correspond to the real solution, is still an issue which needs further studies. 
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